#include "WallGeometry.h"

WallGeometry::WallGeometry(const shared_ptr<CaseSettings> caseSetting1)
{
	//get the K and b for bottom line and top line
	K_bot = ((*caseSetting1).coordinateVertex0[1] - (*caseSetting1).coordinateVertex1[1]) / ((*caseSetting1).coordinateVertex0[0] - (*caseSetting1).coordinateVertex1[0]);
	b_bot = (*caseSetting1).coordinateVertex0[1] - K_bot * (*caseSetting1).coordinateVertex0[0];

	K_top = ((*caseSetting1).coordinateVertex3[1] - (*caseSetting1).coordinateVertex2[1]) / ((*caseSetting1).coordinateVertex3[0] - (*caseSetting1).coordinateVertex2[0]);
	b_top = (*caseSetting1).coordinateVertex3[1] - K_top * (*caseSetting1).coordinateVertex3[0];

	//get normal vector of bottom and top walls
	double middleX = 0.5 * ((*caseSetting1).coordinateVertex0[0] + (*caseSetting1).coordinateVertex1[0]);
	double middleY1= 0.5 * ((*caseSetting1).coordinateVertex0[1] + (*caseSetting1).coordinateVertex3[1]);
	double middleY2 = 0.5 * ((*caseSetting1).coordinateVertex1[1] + (*caseSetting1).coordinateVertex2[1]);
	double middleY = 0.5 * (middleY1 + middleY2);
	vector<double> tempCell{ middleX,middleY};
	normal_bot = GetFaceNorm((*caseSetting1).coordinateVertex0, (*caseSetting1).coordinateVertex1, tempCell);
	normal_top = GetFaceNorm((*caseSetting1).coordinateVertex3, (*caseSetting1).coordinateVertex2, tempCell);
}

double WallGeometry::getBotY(double x0)
{
	return K_bot * x0 + b_bot;
}
double WallGeometry::getTopY(double x0)
{
	return K_top * x0 + b_top;
}

double CalWallDist(const vector<double>& p0, const vector<double>& lineLeft, const vector<double>& lineRight)
{
	// This is 2D method
	// reference: https://math.stackexchange.com/questions/2248617/shortest-distance-between-a-point-and-a-line-segment
	double px, py, vx1, vy1, vx2, vy2;
	px = p0[0]; py = p0[1];
	vx1 = lineLeft[0]; vy1 = lineLeft[1];
	vx2 = lineRight[0]; vy2 = lineRight[1];

	double t1 = -((vx1 - px) * (vx2 - vx1) + (vy1 - py) * (vy2 - vy1)) / (pow(vx2 - vx1, 2.0) + pow(vy2 - vy1, 2.0));
	double dist0 = abs((vx2 - vx1) * (vy1 - py) - (vy2 - vy1) * (vx1 - px)) / sqrt(pow(vx2 - vx1, 2.0) + pow(vy2 - vy1, 2.0));
	double dist1 = min(sqrt(pow(vx2 - px, 2.0) + pow(vy2 - py, 2.0)), sqrt(pow(vx1 - px, 2.0) + pow(vy1 - py, 2.0)));

	if (t1 < 1.0 && t1>0)
	{
		return dist0;
	}
	else
	{
		return dist1;
	}
}

vector<double> GetFaceNorm(vector<double> face_N1, vector<double> face_N2, vector<double> tempCell)
{
	vector<double> tempNormal(3, 0.0);
	face_N1.push_back(0.0);
	face_N2.push_back(0.0);
	tempCell.push_back(0.0);

	vector<double> cell2left{ tempCell[0] - face_N1[0],
	   tempCell[1] - face_N1[1],tempCell[2] - face_N1[2] };

	tempNormal[0] = -1.0 * (face_N2[1] - face_N1[1]);
	tempNormal[1] = face_N2[0] - face_N1[0];

	vector<double> norm2left{ tempNormal[0] - face_N1[0],
	   tempNormal[1] - face_N1[1],tempNormal[2] - face_N1[2] };

	vector<double> right2left{ face_N2[0] - face_N1[0],
	   face_N2[1] - face_N1[1] ,face_N2[2] - face_N1[2] };

	double cellCrossRight = cell2left[0] * right2left[1] - cell2left[1] * right2left[0];
	double normCrossRight = norm2left[0] * right2left[1] - norm2left[1] * right2left[0];
	bool current_direction = true;
	if (cellCrossRight * normCrossRight < 0) // cell center and current norm do not lie on the same side
	{

		tempNormal[1] = -1 * tempNormal[1];
		tempNormal[2] = -1 * tempNormal[2];
	}
	if (cellCrossRight * normCrossRight == 0)
	{
		assert(false && "cell center or normal lies on the face line, not correct!");
	}
	//double dx=localFace.GetNodeID1
	double sqrtNorm = sqrt(tempNormal[0] * tempNormal[0] + tempNormal[1] * tempNormal[1]);
	tempNormal[0] = tempNormal[0] / sqrtNorm;
	tempNormal[1] = tempNormal[1] / sqrtNorm;

	return tempNormal;
}